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1.  INTRODUCTION 

Consistent  and  accurate  coastal  ocean  monitoring  necessitates  the  availability  of 
three  key  components:  (1)  an  observing  network  that  adequately  samples  the  moni¬ 
tored  domain,  (2)  a  coastal  ocean  circulation  model  with  a  sufficiently  high 
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resolution  that  takes  into  account  the  often  complex  geometry  and  dynamics  that 
occur  near  the  coastline,  and  (3)  an  analysis  system  that  is  able  to  accurately  assim¬ 
ilate  the  sampled  observations  to  initialize  the  coastal  model  for  forecasting.  Modern 
analysis  systems  can  also  provide  an  observation  impact  assessment  for  the  design, 
evaluation,  and  possibly  reassignment  of  observing  resources. 

Coastal  current  measurement  types  are  limited  to  moored  buoys  and  ADCPs, 
which  do  not  provide  adequate  spatial  distribution/resolution,  surface  drifters,  which 
tend  to  leave  the  deployment  area  relatively  quickly,  and  shipboard  ADCPs,  which 
are  relatively  expensive  to  operate.  The  continuous  monitoring  of  coastal  waters  for 
circulation  properties  requires  long-term  station  observations.  High-frequency  (HF) 
radar  units  are  unique  observation  platforms  that  provide  surface  current  measure¬ 
ments  at  horizontal  resolutions  of  1  to  6  km  ranging  from  25  to  200  km  off  the  coast. 
This  amount  of  spatial  coverage  would  be  unattainable  with  current  meter  stations. 
In  addition,  HF  radar  units  are  installed  by  various  universities  and  institutions  along 
much  of  the  coastline  of  the  continental  United  States.  HF  radar  observations  have 
been  assimilated  into  ocean  models  mostly  using  sequential  methods  (e.g.,  Refs 
1—4).  However,  there  are  a  few  assimilation  examples  using  a  4DVAR  method.  ~7 

Regional  ocean  models  for  coastal  circulation  monitoring  require  initial  and 
boundary  conditions  from  larger  or  global  domain  models  that  are  usually  run 
with  much  coarser  horizontal  resolution,  as  well  as  surface  forcing  fields  from  atmo¬ 
spheric  models.  To  a  large  extent,  the  accuracy  of  the  coastal  models  depends  on 
(1)  the  accuracy  of  the  larger  domain  model  providing  initial  and  boundary  condi¬ 
tions,  (2)  the  accuracy  of  the  atmospheric  model  providing  surface  forcing  fields, 
and  equally  important,  (3)  the  accuracy  of  the  parameterization  of  the  physics  due 
to  increased  resolution.  The  increased  resolution  can  also  become  a  liability  for 
the  assimilation  as  the  model  resolves  small-scale  circulation  features  that  cannot 
be  constrained  by  the  available  observations.  Usually,  only  coarse  observation 
coverage  is  available  for  assimilation  into  the  larger  domain  (with  the  exception 
of  sea  surface  temperature  (SST)),  making  it  difficult  to  provide  accurate  initial 
and  boundary  conditions  for  the  coastal  model.  Also,  atmospheric  models  can 
contain  errors  in  the  coastal  oceans  due  to  the  coarse  resolution  often  used  and 
the  complex  land— sea  boundary,  not  to  mention  the  lack  of  frequent  feedback 
from  the  ocean  to  the  atmosphere  in  these  areas.  Failure  to  do  this  also  translates 
to  errors  in  the  atmospheric  fields  in  these  areas.  In  addition,  the  ocean  model’s  hor¬ 
izontal  resolution  may  not  be  high  enough  to  capture  all  the  details  of  the  coastline 
and  the  bathymetry.  All  these  elements  contribute  to  the  discrepancies  that  are  seen 
when  coastal  ocean  model  solutions  are  compared  to  observations.  This  is  where  the 
data  assimilation  plays  the  critical  role  of  combining  the  ocean  model  and  available 
observations  in  a  dynamically  consistent  way  to  not  only  provide  a  better  initial  con¬ 
dition  for  the  prediction  of  the  ocean  environment,  but  also  to  correct  at  least  some 
components  of  the  model  error,  e.g.,  errors  in  the  atmospheric  forcing  fields. 

Due  to  the  high  temporal  vaiiability  of  surface  currents  in  the  coastal  areas,  a 
necessary  requirement  of  the  assimilation  system  is  the  ability  to  take  into  account 
the  temporal  dimension  in  the  observations.  Such  capability  is  inherent  to  4DVAR. 
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Contrary  to  the  often  used  sequential  methods  that  assimilate  observations  at  a  given 
time  (thus  correcting  the  model  state  at  fixed  time  stamps),  e.g.,  the  three- 
dimensional  variational  data  assimilation  (3DVAR),  or  the  methods  based  on  the 
Kalman  filter,  4DVAR  seeks  to  correct  the  entire  model  trajectory  for  a  given 
time  window  by  assimilating  all  the  observations  (distributed  in  time  and  space) 
that  were  sampled  during  that  time  window.  In  this  process,  4DVAR  (1)  uses  obser¬ 
vations  at  almost  the  exact  times  that  they  are  sampled,  which  suits  most  asynoptic 
data,  (2)  implicitly  uses  flow-dependent  background  errors,  which  ensures  the  anal¬ 
ysis  quality  for  rapidly  changing  environments,  and  (3)  uses  a  forecast  model  as  a 
constraint,  which  ensures  the  dynamic  balance  of  the  final  analysis. 

It  was  recently  shown  that  the  assimilation  of  surface  velocity  observations 
derived  from  drifters,  using  a  4DVAR  with  the  Naval  Coastal  Ocean  Model 
(NCOM-4DVAR9),  improved  ocean  model  forecasts  of  sea  surface  height,  surface 
and  subsurface  velocity,  temperature,  and  salinity  in  the  Gulf  of  Mexico.9  Unfortu¬ 
nately,  this  study  was  limited  in  time  due  to  the  deployment  and  lifespan  of  the 
drifters.  This  paper  aims  to  expand  on  the  previous  study  by  assimilating  a  sustained 
and  dense  source  of  surface  velocity  observations  from  HF  radar  s  in  the  Chesapeake— 
Delaware  Bay  region  to  show  that  they  are  a  viable  dataset  for  constraining  and 
forecasting  the  coastal  circulation. 


2.  HF  RADAR  OBSERVATIONS 

The  surface  current  observations  used  for  this  study  come  from  a  network  of  three 
SeaSonde  HF  radar  units  that  are  deployed  in  the  mid-Atlantic  region  of  the  East 
Coast  of  the  United  States  (black  stars  on  Figure  1).  The  northernmost  site  is  on 
Assateague  Island,  Maryland,  the  central  site  is  at  Cedar  Island,  Virginia,  and  the 
southernmost  site  is  at  Little  Island  Park,  Virginia.  Throughout  the  study  period, 
data  is  available  from  these  sites  during  93%,  99.9%,  and  100%  of  the  time,  respec¬ 
tively.  These  HF  radar  units  produced  by  CODAR  ocean  sensors  scatter  radio  waves 
off  the  ocean  surface  and  infer  movement  of  near  surface  currents.  During  July  20 1 3, 
the  three  stations  were  operating  at  4.5  MHz  that  resonantly  scatter  off  surface  grav¬ 
ity  waves  of  approximately  30-m  wavelength.  For  these  so-called  long-range  site 
observations  are  provided  hourly  on  a  polar-coordinate  grid  with  a  range  step  of 
6  km  and  a  bearing  step  of  5°.  Additionally,  the  horizontal  range  of  a  single  site 
is  approximately  200  km  offshore. 

Roarty  et  al.10  discusses  the  operation  and  maintenance  of  the  mid-Atlantic  HF 
radar  observing  network;  this  includes  the  three  sites  used  in  this  study.  Because  the 
HF  radar  data  is  used  by  the  U.S.  Coast  Guard  Search  and  Rescue  Optimal  Planning 
System  (S  AROPS),  there  is  an  implemented  procedure  for  quality  control  and  assur¬ 
ance  of  the  observations  collected  by  these  sites.  Additionally,  they  go  on  to  report 
RMS  differences  between  in  situ  measurements  of  both  acoustic  doppler  current  pro¬ 
filers  (ADCPs)  and  drifters  of  7.4  to  9.8  cm/s  using  an  unweighted  least  squares 
method. 1 1  This  unweighted  least  squares  approach  merges  the  single-site  radial 
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The  model  domain:  the  outer  box  is  the  3-km  parent  nest,  and  the  inner  and  expanded  box  is 
the  actual  1-km  domain  with  a  sample  coverage  of  HR  radar  observations  on  July  01,  2013, 
overlaid  on  the  colored  bathymetry.  The  radar  stations  are  represented  by  the  stars. 

velocities  located  within  a  search  radius  around  each  grid  point.  This  processing  step 
uses  a  Matlab  toolbox  called  ITFRJProgs  to  create  total  current  vectors.  For  the 
three-site  setup  used  here,  the  search  radius  is  10  km  with  a  minimum  of  two  radials 
from  at  least  two  sites  required  to  create  a  total.  In  theory,  accurate  surface  velocities 
are  recovered  when  two  HR  radar  beams  form  an  angle  of  90°.  In  practice,  however, 
accurate  velocities  can  still  be  constructed  from  beams  forming  an  angle  as  low  as 
15°,  which  has  become  a  standard  for  operational  processing  of  HF  radar  observa¬ 
tions,  at  least  in  the  mid-Atlantic  Bight  HF  radar  observing  network. 111  The  HF  radar 
observations  assimilated  here  were  processed  with  the  15°  minimum  angle 
threshold.  The  HF  radar  measurements  are  sensitive  to  environmental  factors  that 
can  affect  the  spatial  extent  of  the  velocity  footprint.  This  usually  results  in  occa¬ 
sional  gaps  within  a  coverage  area.  For  more  information  see  Ref.  12. 


3.  THE  MODEL 

The  ocean  model  used  in  this  study  is  the  Navy  coastal  ocean  model  (NCOM). 
NCOM  is  a  free-surface  model  that  has  been  described  in  the  literature.1  14  The 
model  domain  (and  bathymetry)  shown  in  Figure  1  spans  longitudes  76.6°W  to 
73. 6° W  and  latitudes  36.75°N  to  39.2°N  at  1-km  horizontal  resolution  with  50 
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vertical  levels.  The  initial  conditions  were  obtained  from  downscaling  the  opera¬ 
tional  1/8°  resolution  global  NCOM  (GNCOM)  to  an  intermediate  model  with  hor¬ 
izontal  resolution  of  3  km,  and  then  to  a  high-resolution  1  -km  model.  Horizontal 
viscosities  and  diffusivities  are  computed  using  either  the  grid-cell  Reynolds  number 
(Re)  or  the  Smagorinsky  schemes,  both  of  which  tend  to  decrease  as  resolution  is 
increased.  The  grid-cell  Re  scheme  sets  the  mixing  coefficient  K  to  maintain  a 
grid  cell  Re  number  below  a  specified  value,  e.g.,  if  Re  =  u  *  dx/K  =  30,  then 
K  =  u  *  dx/30.  Hence,  as  dx  decreases,  K  decreases  proportionally.  A  similar 
computation  is  performed  for  the  Smagorinsky  scheme. 

The  surface  atmospheric  forcing,  including  wind  stress,  atmospheric  pressure, 
and  surface  heat  flux,  is  provided  by  the  Navy  Global  Atmospheric  Prediction  Sys¬ 
tem  (NOGAPS l:’-17)  with  a  horizontal  resolution  of  0.5°.  River  forcing  is  provided 
at  all  river  in-flow  locations  in  this  mid-Atlantic  domain.  Additionally,  eight  tidal 
constituents  (Kl,  01,  PI,  Ql,  K2,  M2,  N2,  and  S2)  are  forcing  the  domain  through 
the  open  boundaries.  Open  boundary  conditions  use  a  combination  of  radiative 
models  and  prescribed  values  provided  by  the  parent  3-km  nest.  Different  radiative 
options  are  used  at  the  open  boundaries  depending  on  the  model  state  variables:  a 
modified  Orlanski  radiative  model  is  used  for  the  tracer  fields  (temperature  and 
salinity),  an  advective  model  for  the  zonal  velocity  (u),  a  zero  gradient  condition 
for  the  meridional  velocity  (v)  as  well  as  the  barotropic  velocities,  and  the  Flather 
boundary  condition  for  elevation. 


4.  THE  ASSIMILATION  SYSTEM 

The  assimilation  system  used  here  is  described  in  more  detail  in  Ref.  8.  The  brief 
presentation  that  follows  only  serves  to  elucidate  the  focus  of  this  study.  For  a  given 
model,  the  following  is  presented: 

|f  =  F(x)+/,o<,<r  (1) 

[  x(t  =  0)  =  /(*)  +  i(x) 

where  X  stands  for  all  the  dependent  model  state  variables,  i.e.,  the  two-dimensional 
SSH  and  barotropic  velocities,  and  the  three-dimensional  temperature,  salinity,  and 
baroclinic  velocities;  F  includes  the  model  tendency  and  forcing  terms,  /  is  the 
model  error  with  covariance  Cf,  I(x)  is  the  prior  initial  condition,  and  i{x)  is  the  initial 
condition  error  with  covariance  C,;  x  and  t  represent  the  position  in  the  three- 
dimensional  space  and  time,  respectively.  Given  a  vector  Y  of  M  observations  of 
the  model  state  in  the  space— time  domain,  with  the  associated  vector  of  observation 
errors  e  (with  covariance  CE),  the  following  is  shown: 


y,n  —  HmX  +  em 


1  <  m  <  M 


(2) 
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where  Hm  is  the  observation  operator  associated  with  the  wth  observation.  One  can 
define  a  weighted  cost  function  as  follows: 


T  T 

J  =  j  II  j  f{x,  t)Wf(x,  t,x,  t')f(x  ,t')dx'dt'dxdt 
o  b  o  a 


+ 


I  I  i(x)Wt(x,x')i(x!)dx 

h  o 


'dx  +  eTWee 


(3) 


where  Q  denotes  the  model  domain,  the  weights  Wf  and  W,  are  defined  as  inverses  of 
Cf  and  C,  in  a  convolution  sense,  and  Wt  is  the  matrix  inverse  of  Cc.  Boundary 
condition  errors  are  omitted  from  Eqns  ( 1 )  and  (3)  only  for  the  sake  of  clarity. 
The  solution  of  the  assimilation  problem,  i.e.,  the  minimization  of  the  cost  function 
(Eqn  (3)),  is  achieved  by  solving  the  following  Euler-Lagrange  (EL)  system: 


—  =  F(X)  +  Cr X,  0<t<T, 

X(t  =  0)  =  I{x)  +  CfX{xt0) 

T  MM 

^  W^mn  (ym  ~  H"’x)d(x  -  xm)6(t  -tm),  0  <  t  <  T 


dX 

~di 


d> 


m= 1  n— 1 


A(r)  =  o 

where  X  is  the  adjoint  variable  defined  as  the  weighted  residual: 


(4) 


T 

X(x,t)  =  J  I  Wf(x,t,x',t')f{x',t,)dx'dt',  (5) 

o  b 

and  <5  denotes  the  Dirac  delta  function,  Wf:nm  are  the  matrix  elements  of  the 
superscript  T  denotes  the  transposition,  and  the  covariance  multiplication  with  the 
adjoint  variable  is  the  convolution: 


and 


r 

CfX(x,t)=  I  I  Cf(x,t,x',t')X(x/,t')dx'dt', 
0  b 


C;°A(.t,0)=  I  Ci(x,x)X(x'.Q)dx' 

b 


(6) 


(7) 


for  the  model  and  initial  condition  errors,  respectively. 
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It  can  be  seen  in  Eqn  (4)  that  the  adjoint  model  is  forced  by  the  innovations 
(model-data  misfits  at  the  observation  locations),  and  its  solution  initializes  and/or 
forces  the  forward  model,  depending  on  whether  a  strong  or  weak  constraints 
assumption  is  adopted. 

A  standard  approach  to  solving  the  Euler-Lagrange  system  (Eqn  (4))  is  the 
strong  constraints  4dvar  that  assumes  that  only  the  initial  condition  is  erroneous, 
i.e.,  the  model  has  no  errors  (C/  =  0).  The  solution  of  Eqn  (4)  is  found  iteratively 
as  follows:  (1)  a  first  guess  initial  condition  is  used  to  solve  the  nonlinear 
model,  (2)  the  nonlinear  solution  is  used  to  compute  the  model-data  misfits  that 
appear  in  the  right-hand  side  of  the  adjoint  model,  (3)  the  adjoint  model  is  solved 
and  used  to  compute  the  correction  to  the  initial  condition,  and  (4)  the  process  is 
repeated  until  the  minimum  of  the  cost  function  or  a  preselected  convergence  cri¬ 
terion  is  reached. 

The  weak  constraints  4dvar  approach  takes  into  account  the  model  errors  and, 
thus,  increases  the  dimension  of  the  control  space,  which  now  becomes  the  entire 
model  trajectory  for  the  selected  assimilation  window.  This  rather  huge  control 
space  also  increases  the  computational  cost  of  the  assimilation,  and  it  usually  renders 
the  minimization  (of  the  cost  function)  process  poorly  conditioned.  This  difficulty 
can  be  avoided  if  the  minimization  is  done  in  the  data  space,  which  does  not  depend 
on  and  is  usually  much  smaller  than  the  control  space.  That  is  possible  through  the 
representer  algorithm, 1 8' 19  which  expresses  the  solution  of  Eqn  (4)  as  the  sum  of  a 
first  guess  and  a  finite  linear  combination  of  representer  functions,  one  per  datum. 
Being  a  linear  expansion,  the  representer  algorithm  cannot  be  applied  to  Eqn  (4) 
directly,  mainly  because  of  its  nonlinear  property.  However,  following  Refs  8,20 
the  representer  algorithm  can  be  applied  to  a  linearized  form  of  Eqn  (4). 


5.  EXPERIMENTS  AND  RESULTS 

For  this  study,  the  initial  condition  error  for  the  experiment  is  set  as  1.0  °C  for  tem¬ 
perature,  0.1  practical  salinity  unit  (psu)  for  salinity,  and  0.5  m/s  for  velocity.  These 
errors  are  set  by  examining  the  innovation  values  between  a  free-running  NCOM 
and  available  observations.  The  error  values  are  uniformly  prescribed  across  the 
domain  and  reduced  at  depth.  This  is  deemed  acceptable  because  we  are  mostly 
interested  in  the  accuracy  of  surface  currents.  It  is  also  important  to  note  the  model 
errors  in  this  study  are  attributed  to  errors  in  the  specified  atmospheric  surface  forc¬ 
ing.  This  a  reasonable  assumption  as  ocean  surface  currents  are  strongly  influence  by 
surface  wind  stress.  The  model  errors  are  0.05  °C  for  temperature,  0.005  practical 
salinity  unit  (psu)  for  salinity,  and  0.05  m/s  for  velocity.  The  model  errors  represent 
5%  of  the  magnitudes  of  the  atmospheric  forcing  in  respective  equations,  with 
the  exception  of  the  free  surface,  and  are  converted  from  fluxes  units  to  units  of 
the  ocean  state  variables  using  the  relationships  imposed  by  the  discretization 
of  the  model  (see  Refs  21,8).  The  horizontal  correlation  scales  of  the  initial  and 
model  errors  are  taken  to  be  20  km  (approximately  the  Rossby  radius  of  deformation 
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in  the  domain)  and  fixed  in  time.  When  these  isotropic  covariances  are  convolved 
with  the  adjoint  solution  and  included  in  the  forward  model  as  dictated  by  Eqn 
(4),  the  4dvar  system  produces  analysis  increments  that  are  flow  dependent,  thanks 
to  the  dynamics  of  the  tangent  linear  and  adjoint  models,  and  also  the  nonlinear 
model  trajectory  around  which  the  system  is  linearized. 

Each  of  the  observation  data  types  is  also  assigned  errors.  The  observation  errors 
are  usually  a  combination  of  the  estimated  instrument  error  and  the  representative¬ 
ness  error.  Here,  the  temperature  error  is  0.35  °C,  0.035  practical  salinity  unit  (psu), 
and  0.05  m/s  for  velocity.  The  experiment  carried  out  here  takes  place  from  July  1  to 
31, 2013,  in  sequential  assimilation  windows  of  three  days,  with  observations  binned 
hourly  and  subsampled  to  keep  only  one  observation  per  20-km  correlation  scale  in 
both  meridional  and  zonal  directions.  With  the  exception  of  the  first  cycle,  the  back¬ 
ground  (i.e.,  the  solution  that  the  assimilation  is  trying  to  correct)  for  each  cycle  is 
the  forecast  obtained  by  running  the  nonlinear  with  the  final  condition  from  the  anal¬ 
ysis  in  the  previous  cycle. 

In  order  to  assess  the  fit  to  the  observations  over  time  in  the  whole  assimilation 
window,  we  define  the  following  “fit  to  the  observations”  metric: 


(8) 


In  Eqn  (8),  ym  is  the  with  observation,  M  is  the  total  number  of  observations,  H,„  is 
the  observation  operator,  X“  is  the  assimilated  solution  or  analysis,  and  am  is  the 
observation  error  or  standard  deviation.  The  right-hand  side  of  Eqn  (8)  can  be 
computed  as  a  time  series  and  also  evaluated  for  the  free-run  solution  and  the  first 
guess.  Because  the  assimilation  is  expected  to  fit  the  observations  to  within  the 
observation  standard  deviation  at  the  observation  locations,  the  metric  Jpij-  in 
Eqn  (8)  is  expected  to  be  less  or  equal  to  one  for  the  analysis.  One  only  hopes 
that  the  same  is  true  for  the  subsequent  forecasts  as  a  result  of  fitting  the  observations 
in  previous  cycles. 

The  results  of  this  assimilation  experiment  show  that  the  NCOM-4DVAR  is 
capable  of  assimilating  HF  radar  velocities  by  significantly  reducing  the  discrep¬ 
ancies  between  the  modeled  and  the  observed  surface  velocities.  It  can  be  seen 
in  Figure  2  that  the  free-running  model  for  the  1  -km  resolution  (black  line)  is  in  sig¬ 
nificant  disagreement  with  the  observations,  having  Jpn  values  between  2  and  4 
observation  standard  deviations.  The  assimilation  (gray  dashed  line)  is  able  to 
reduce  those  discrepancies,  sometimes  by  as  much  as  2  standard  deviations,  with 
Jfir  values  generally  between  1  and  2.4  observation  standard  deviations.  These 
values  are  still  higher  than  the  target  value  of  1 ,  indicating  that,  although  the  assim¬ 
ilation  has  done  a  good  job  of  reducing  the  discrepancies  compared  to  the  free- 
running  model,  the  assimilated  solution  is  still  not  fitting  the  observations  accurately. 
On  the  other  hand,  the  first  guess  solution  (dashed  black  line),  which  consists  of  the 
forecast  from  analysis  in  the  previous  assimilation  window,  shows  discrepancies  of 
the  same  magnitudes  as  the  free-running  solution.  This  indicates  that  the  gains  of  the 
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FIGURE  2 


The  Jfit  metric  for  a  free-running  model  (black  dashed  line),  the  4DVAR  first  guess  fields 
(solid  black  line),  and  the  4DVAR  analysis  (gray  dashed  line)  at  the  observation  locations. 
This  experiment  assimilates  observations  hourly. 


assimilation  are  quickly  lost  in  the  forecast,  primarily  due  to  erroneous  surface  forc¬ 
ing  and  the  high  resolution  of  the  model  that  resolves  circulation  features  that  are  not 
observed. 


6.  VALIDATION 

The  validation  of  any  assimilation  experiment  requires  independent  observations 
against  which  the  assimilation  results  can  be  compared.  Those  usually  consist  of 
buoys  along  the  coast.  However,  those  are  not  available  during  the  time  of  this  exper¬ 
iment.  We  carried  out  a  second  assimilation  experiment  where  observations  are 
assimilated  every  3  h  instead  of  every  hour  as  in  the  previous  experiment.  The  un¬ 
assimilated  observations,  i.e.,  those  that  are  left  out  every  2  h,  are  considered  inde¬ 
pendent  for  the  purpose  of  validation,  even  though  they  may  be  correlated  with  those 
that  are  assimilated,  by  reason  of  proximity  in  space  and  time.  The  same  fit  to  the 
observations  metric  J in  is  also  used  to  evaluate  this  assimilation  experiment,  not 
only  for  the  assimilated  observations,  but  also  for  the  withheld  observations.  Results 
in  Figure  3(a)  show  the  JF/T  values  for  the  assimilated  solution  and  the  first  guess 
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FIGURE  3 

Jor  metric  for  the  first  guess  solution  (dashed)  and  the  new  assimilated  solution  (solid)  at  the 
assimilated  observations  every  third  hour  (a)  and  at  the  withheld  observations  (b). 
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compared  to  assimilated  observations,  whereas  Figure  3(b)  shows  similar  .////  values 
for  the  assimilated  solution  and  the  first  guess  compared  to  withheld  observations.  It 
can  be  seen  that  similar  to  the  previous  experiment,  there  is  a  significant  reduction  in 
the  Jfir  values  from  the  first  guess  (2—3.5  standard  deviations)  to  the  assimilated 
solution  (1.2— 2.7  standard  deviations).  More  importantly,  there  is  a  good  improve¬ 
ment  of  the  assimilated  solution  versus  the  first  guess  when  these  two  solutions  are 
compared  to  the  withheld  observations,  an  improvement  that  sometimes  exceeds  a 
standard  deviation. 

Figure  4  shows  a  comparison  of  surface  velocities  maps  from  the  observations, 
the  first  guess,  and  the  analysis  5  and  16  days  into  the  assimilation,  respectively. 
We  first  note  that  at  these  two  time  levels,  there  is  almost  no  agreement  at  all 
between  the  circulation  patterns  shown  in  the  observations  and  those  in  the  first 
guess,  i.e.,  the  model  is  significantly  in  error  compared  to  the  observations,  even 
after  being  initialized  by  the  assimilation.  For  example,  on  day  5,  the  observations 
describe  an  offshore  surface  circulation,  whereas  the  model  shows  an  alongshore 
circulation.  This  results  from  the  atmospheric  forcing  fields  being  erroneous 
themselves,  see  Figure  5.  The  assimilation  procedure  alters  the  first  guess  enough 
to  produce  an  analysis  that  fits  (looks  like)  the  observations,  albeit  not  perfectly: 
The  offshore  current  is  reconstructed  by  the  assimilation  on  day  5;  and  on  day  16, 
the  northeastward  coastal  current  that  turns  offshore  is  also  recovered,  though  these 
features  were  missing  in  the  first  guess.  However,  the  analysis  sometimes  still  has 
the  patterns  of  the  first  guess.  This  indicates  that  the  background  and  model  errors 
prescribed  for  the  assimilation  are  too  small. 

A  major  difficulty  in  coastal  ocean  modeling  resides  in  the  lack  of  high  spatial 
resolution  atmospheric  forcing.  Atmospheric  models  are  usually  run  with  coarser 
resolutions  compared  to  the  ocean  models  (especially  for  coastal  applications), 
because  resolving  the  rather  fast  motion  of  the  atmosphere  with  high  horizontal 
resolution  would  require  very  small  time  steps  that  would  be  computationally 
prohibitive.  For  the  case  at  hand,  the  ocean  model  has  a  horizontal  resolution 
of  1  km,  while  the  atmospheric  forcing  fields  are  obtained  from  interpolating 
results  from  an  atmospheric  model  that  used  a  0.5°  resolution  that  does  not  capture 
the  variability  of  the  model  domain.  According  to  Ekman  theory,  a  modest  wind 
stress  of  5  m/s  would  cause  the  surface  velocity  to  deflect  to  the  right  of  the 
wind  stress  direction  by  an  angle  of  45°  if  the  water  depth  is  at  least  45  m.  The 
topography  shown  in  Figure  1  indicates  that  this  would  not  apply  to  a  significant 
portion  of  the  domain  where  the  depth  is  less  than  45  m,  and  it  is  expected  that 
the  surface  velocity  be  strongly  correlated  to  the  wind  stress  in  that  part  of  the 
domain.  Figure  5  shows  the  direction  of  the  wind  stress  compared  to  the  direction 
of  the  surface  currents  from  the  HF  radar  stations.  The  wind  stress  is  generally  uni¬ 
form  as  a  result  of  interpolating  from  a  few  grid  points  of  the  atmospheric  model.  It 
can  be  seen  that  the  direction  of  the  wind  is  never  aligned  with  the  direction  of  the 
observations;  they  are  quite  different.  This  presents  a  significant  challenge  to  the 
assimilation  system  and  explains  why  the  assimilation  could  not  accurately  fit 
the  observations  with  small  initial  conditions  and  model  errors. 


Surface  velocities  from  the  observations  (a,  d),  the  first  guess  (b,  e),  and  the  analysis  (c,  f)  5  and  16  days  into  the  assimilation,  respectively. 
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A  comparison  of  the  direction  of  the  surface  currents  (a,  c,  and  e)  and  the  direction  of  the 
wind  stress  (b,  d,  and  f)  at  the  end  of  the  first,  the  sixth,  and  the  tenth  assimilation  window, 
respectively. 


Another  challenge  to  the  assimilation  resides  in  the  resolution  of  the  ocean  model 
and  that  of  the  observations.  The  processed  observations  have  a  spatial  resolution  of 
about  6  km,  whereas  the  model  is  run  at  a  resolution  of  1  km.  As  seen  in  Figure  6  in  a 
subset  of  the  domain,  the  model  resolves  multiple  small-scale  features  that  are  ab¬ 
sent  from  the  observations.  Also,  the  strong  flow  at  the  southeast  corner  of  the 
domain  reveals  that  the  model  boundary  conditions  contain  an  intrusion  of  the 
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FIGURE  6 


A  comparison  of  the  1-km  first  guess  velocity  field  (a)  and  the  observations  (b)  on  day  10  in  a 
subset  of  the  model  domain. 


Gulf  Stream,  another  feature  that  is  not  present  in  the  observations.  Thus,  the  assim¬ 
ilation  could  benefit  from  more  accurate  boundary  conditions  and,  perhaps,  a 
slightly  coarser  resolution  from  the  model,  say  3  km. 

A  third  assimilation  experiment  is  carried  out  for  three  cycles  of  3  days  each, 
with  the  same  setup  as  the  original  experiment,  except  that  the  model  resolution 
is  reduced  from  1  to  3  km.  Figure  7(a)  shows  a  comparison  between  the  analyses 
of  the  1  and  3  km  experiments,  where  significant  improvements  in  the  accuracy  of 
the  3-km  analysis  can  be  seen,  especially  at  the  times  when  the  1-km  analysis  has 
J fit  values  exceeding  1 .5  standard  deviations.  In  general,  JFn  values  for  the  3-km 
analysis  are  lower  than  1 .5  standard  deviations.  On  the  other  hand,  a  similar  compar¬ 
ison  of  J/.yy  values  for  the  first  guesses  from  both  1-  and  3-km  experiments  in 
Figure  7(b)  shows  that  the  forecast  from  the  3-km  analysis  has  significantly 
improved  compared  to  the  1-km  forecast,  with  noticeable  gaps  between  the  two  lines 
sometimes  exceeding  1  standard  deviation.  However,  the  3-km  first  guess  still  dis¬ 
plays  large  discrepancies  with  the  observation,  having  Jpn  values  generally 
exceeding  2  observation  standard  deviations,  and  only  occasionally  falling  below 
1.5  standard  deviations.  Once  again,  this  loss  of  accuracy  in  the  first  guess  is  attrib¬ 
uted  to  the  erroneous  atmospheric  forcing  (in  this  case  the  wind  stress)  because  the 
model  resolution  is  now  closer  to  that  of  the  observations  coverage,  and  the  analysis 
is  also  significantly  closer  to  the  observations. 

Similar  to  Figure  4,  Figure  8  shows  a  comparison  of  surface  velocities  maps  from 
the  observations,  the  first  guess,  and  the  analysis  3  and  5  days  into  the  assimilation, 
respectively,  for  the  3-km  model  resolution.  There  is  a  significant  difference  in  the 
direction  of  the  velocity  between  the  observations  and  the  first  guess  on  July  3,  espe¬ 
cially  at  the  lower-right  side  of  the  domain  showing  an  intrusion  of  the  Gulf  Stream 
in  the  first  guess,  whereas  the  observations  locate  this  intrusion  further  to  the  east 
and  slightly  to  the  north,  compared  to  its  location  in  the  first  guess.  The  assimilation 
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FIGURE  7 


A  comparison  of  the  JF,T  values  from  the  analyses  (a)  and  the  first  guesses  (b)  from  the  1-km 
(dashed)  and  3-km  (solid)  model  resolutions. 


Similar  to  Figure  4,  surface  velocities  from  the  observations  (a,  d),  the  first  guess  (b,  e),  and  the  analysis  (c,  f)  on  July  3  and  July  5,  respectively. 
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procedure  corrects  the  first  guess  significantly  to  produce  an  analysis  that  fits  (looks 
like)  the  observations,  e.g.,  the  Gulf  Stream  current  is  in  better  agreement  with  the 
observations,  albeit  not  perfectly.  However,  the  analysis  on  July  3  still  has  some  pat¬ 
terns  of  the  first  guess,  indicating  that  the  background  and  model  errors  prescribed 
for  the  assimilation  may  still  be  too  small.  The  analysis  on  July  5  shows  a  much  bet¬ 
ter  agreement  with  the  observations  (offshore  circulation),  even  though  the  first 
guess  was  not  (coastal  circulation).  Thus,  even  in  the  presence  of  an  erroneous 
wind  stress,  the  weak  constraint  assimilation  fits  the  observations  significantly  better 
when  the  model’s  horizontal  resolution  is  closer  to  that  of  the  observations  coverage. 


7.  CONCLUSION 

Surface  velocity  observations  from  HF  radar  are  a  valuable  data  set  for  monitoring 
the  coastal  circulation,  and  they  can  be  assimilated  into  a  coastal  circulation  ocean 
model  using  the  NCOM-4DVAR  system,  provided  adequate  model  resolution, 
initialization,  boundary  conditions,  and  atmospheric  forcing.  It  was  shown  in  the 
experiments  presented  here  that  the  assimilation  cannot  accurately  fit  the  observa¬ 
tions  with  rather  small  initial  conditions  and  model  errors.  However,  the  biggest 
challenge  for  the  assimilation  system  consists  in  an  erroneous  wind  stress  that 
consistently  steers  the  model  in  a  completely  different  direction  than  the  observed 
surface  velocities.  Although  the  assimilation  was  able  to  reduce  a  noticeable  portion 
of  the  model’s  discrepancy  to  tire  observations,  those  gains  were  quickly  lost  in  the 
forecast  stage  for  the  following  assimilation  window,  primarily  due  to  the  wrong 
wind  stress  and  the  high  resolution  of  the  model.  Reducing  the  model  resolution 
to  be  closer  to  that  of  the  observations  significantly  improved  the  accuracy  of  the 
analysis  and  the  forecast.  The  ability  of  the  assimilation  system  to  accurately  fit 
the  observations  can  also  be  improved  by  prescribing  higher  model  errors.  However, 
we  suggest  that  instead  of  increasing  the  error  levels  in  the  assimilation  system, 
which  can  be  justified  for  the  case  at  hand,  the  primary  source  of  the  errors  must 
be  addressed  first,  i.e.,  providing  a  wind  stress  that  drives  the  model  to  be  in 
more  agreement  with  the  observations.  This  can  be  achieved  by  a  local  nest  of  the 
atmospheric  model,  or  better  yet,  a  fully  coupled  ocean— atmosphere  model. 
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